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Abstract 

This interim report describes the research being conducted in the formulation of 
hierarchic models for laminated plates. The work, conducted in collaboration with 
Professors Babuska and Schwab of the University of Maryland, is an extension of 
the work done for laminated strips [3] . The use of a single parameter /?, representing 
the degree to which the equilibrium equations of three-dimensional elasticity are 
satisfied, is being investigated. The powers of f3 identify members in the hierarchic 
sequence. 

Included in this report are numerical examples that were analyzed with the 
proposed sequence of models. The results obtained for square plates with uniform 
loading and with homogeneous boundary conditions are very encouraging. Several 
cross-ply and angle-ply laminates were evaluated and the results compared with 
those of the fully three-dimensional model, computed using MSC/PROBE, and 
with previously reported work on laminated strips [3]. 
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1 Models for laminated plates 


Consider an infinite flat plate of constant thickness h composed of thin layers of 
orthotropic material perfectly bonded together. Each layer (lamina) possesses a 
plane of elastic symmetry parallel to the x—y plane. The laminae are symmetrically 
arranged with respect to the middle surface of the plate (i.e., the x — y plane). The 
load q(x, y) is antisymmetric with respect to the middle plane, and q(x, y) = 0 for 
M > a, \y\ > 6, with a and 6 some fixed number. Let a = 1 /a and 7 = 1/6, and 
further let: 

P = min(a, 7) (1) 

Assume that the displacement field can be written in exponential form: 

u x (x,y,z) = 4 >{( 3 , z) e ,/J(r+3,) (2) 

u y {x,y,z) = 0 (/ 3 , A e * P( x +v) ( 3 ) 

u z {x,y,z) = p(ft, z) e'^ x+y) ( 4 ) 

where 

0 (/?,z) = 4 > a (/ 3 ,z)+i<j) b {/ 3 ,z) ( 5 ) 

0(A z) = 0a (A z)+iifo(P, z) (6) 

p(P, z) - pa(fi,z)+ipb{ft, z) ( 7 ) 

where 0 a , 0„, 06 are antisymmetric real functions, and p a and p b are symmetric 

real functions with respect to the middle surface of the plate (laminate) . 

The strain components corresponding to the displacement field given by (2), ( 3 ) 
and ( 4 ) are: 
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where the primes represent differentiation with respect to z. 

Let x\ y' be the material (lamina) coordinates for any layer rotated an angle 
6 with respect to the global coordinate system about the z axis. Then the stress- 
strain relations in the global (laminate) system for any layer can be written as: 




(14) 


where [C] is the lamina material stiffness matrix in the lamina coordinate system 
(x\ y', z), and [T] is the transformation matrix. Defining 

[Q] = [T\~ 1 [C\[T\ (15) 

as the transformed lamina material matrix, equation (14) can be written as: 
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The equilibrium equations with zero body force components are given by: 

dcr x 


dx 
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dr yz dcr z _ „ 
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(19) 


Substituting equations (8) to (13) and (16) into (17) to (19) the following equi- 
librium equations are obtained: 


{ — /? 2 [(<5ll + 2Qi6 + Qte)<}> + (Qn + Ql6 + Q 26 + Q66)'0] + 
* /?[(Qi3 + Qi^p' + (Q 45 p y + (Q 55 p)'\ + 
(Q45^ +Qss<f>'Y}e i ^ x+9) =0 


(20) 
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{ — /? 2 [(Ql2 + Ql 6 + Q26 + + (Q22 + 2Q26 + + 

i fl[{Q 23 + Q 3 g)p' 4“ (Q44 py + (Q45 p) 1 ] + 

(<?44V' , + g4 5 0 / ) , }e ,/?(l+J ' ) =O (21) 

{—(^{{Qa + 2 Q 45 + Qhh)p\ + 

* /3[(Q44 + QuW + (Q 45 + Qtt)# + ((Ql3 + QnW + 

((Q 23 + QzsWl 4- (Q33 p')'}e'^ xJry ^ = 0 (22) 

Expanding 0(/3, z), z) and p((3, z) into a Taylor series with respect to (3: 

(f>((3, z) = [(f) a 0 (z) + i (pbo (z)\ 4- j3[(f>ai (z) + i 4> bx (z)} + (3 2 [<j> a2 (z) + i (f) b2 {z)] 4 

1p(/3, z) = [i>ao(z) +ilko(z)] +Pbl>al{z) + i^l(z)] + fi 2 [Mz) +* W*)] +** 1 
p{/3, z) - [pao{z)+ipbo{z)\+/3[pal{z)+ipbl{z)\+(3 2 [pa2{z)+ipb2{z)\-\ 

On substituting into the equilibrium equations (20), (21), (22) and separating 
into real and imaginary parts (only the real parts are shown here) we have: 

The real part of (20): 

{Q^aO + Q^4>ao)' + AHQiSTpal + Qttfial)' ~ (Ql3 + Q^PbO ~ 

(( Q45 + Q5s)pbo)'] + 

+ + Qttfial)' ~ {Q\Z + Q 26 )Pbl ~ 

((^45 + Qb^Pbl)' ~ ( Qll +2Qi6 + ^66)^0 — 

(Ql2 + Ql 6 + Q26 + ^ 66 )^ 00 ] 4 =0 (23) 

The real part of (21): 

{Qu^aO + Qte^ao)' + /5[(Q440ai + Qttftal)' ~ (Q» + QzsjPbQ ~ 

(($44 + Q4s)/26o) T + 

+ /^ 2 [(Q44'0a2 + Qttfirt)' ~ ($23 + Q36)Pbl ~ 

(($44 + Q4s)Pbl)' ~ (Ql2 + $16 + $26 + $66)<&i0 ~ 

(Q22 + 2Q26 + $66)^0] 4 — 0 (24) 

The real part of (22): 

(Q 33 PaoY + ft[(Q 33 pal)' ~ (Q44 + $45)V , &> “ ($45 + QttWbO ~ 

(($13 + Q36)<f>bo)' ~ (($23 + Q36)Vfa)) / ] + 

4* (P[{Q 33 Pa 2 )' ~ ($44 + $ 45)^1 ~ ($45 + QttWbl ~ 

(($13 4" Q3&) < f^b\) > ~ (($23 +Q3f)' l k>\)' ~ 

($44 4- 2 Q 45 4- $ 5 s)Aio] 4 =0 (25) 

These equations hold for any choice of f3. Solving for each power of [3 we obtain 
the transverse shape functions as described in next section. 


3 


1.1 The model characterized by /3° 

Setting (3 = 0 in equations (23) to (25) we have: 


( Q^'aO + Q^ao)' — 0 (26) 

(Qu^+QM = 0 (27) 

(QkP’oY = 0. (28) 


Knowing that 4> a 0 (z), ip a o(z) are antisymmetric and p a 0 (z) is symmetric, and inte- 
grating, we have: 


M*) = Fo(z) (29) 

Tpao{z) = G 0 (z) (30) 

pao(z ) = 1. (31) 


Similarly, solving the imaginary part of the equilibrium equations we get: 


4 >bo(z) = F 0 {z) 
i>bo{z) = Go(z) 

Pbo(z ) = 1 . 


where 


F a {z) 


G 0 (z) 


r Qu — Q45 
/ QuQss - Q45 

f Qss ~ Q45 
J Q44Q55 - QI5 


dz 


dz. 


(32) 

(33) 

(34) 


(35) 

(36) 


To obtain these transverse functions we adopted all integration constants to be 
either 0 or 1. This accomplishes one important aspect, i.e., there is only ‘one’ 
transverse function per field. If the integration constants are arbitrary, the solution 
is: 


<t>ao(z) 


lpao(z) 


r op Q44 — bp Q45 

J QuQss — Qh 


dz 


r bp Q55 — a o Q45 
J Q44Q55 - Qh 


dz 


(37) 


(38) 



Note that with either choice of the integration constants, the resulting transverse 
functions (29), (30), (31) or (37), (38), (39) satisfy the equilibrium equations. Se- 
lecting the integration constants to be some convenient numbers is acceptable be- 
cause the equilibrium equations are still satisfied. 

The real and imaginary parts are not linearly independent, hence both lead to 
the same functional form The mode of deformation corresponding to /? = 0 can 
be written in the following form 

u x (x,y,z) 
u y (x,y,z) 
u z {x,y,z) 

When $ 44 , Q45, $55 are constant trough the thickness, this model is capable of 
representing rigid body displacement and rotation. However, this model does not 
satisfy the condition of converging to the same limit as the problem of elasticity 
as h — *• 0, unless some adjustments are introduced to the material properties as 
discussed later. 


= Ui(x, y)Fa(z) 

(40) 

= u 2 {x, y)Go(z) 

(41) 

= u 3 (x,y). 

(42) 


1.2 The model characterized by j3 l 

To find the mode of deformation for the model which satisfies the equilibrium 
equations up to the first power of /?, we differentiate (23) to (25) with respect to j3 
and let /? = 0. In this case we have: 

($450ol + $55001 y ~ ($13 + Q^PbO ~ (($45 + $ 55 )^ 60 ) 7 := 0 (43) 

($440al 4* $450oi) 7 — ($23 + $36)^60 — (($44 + $ 45 )^ 60 ) ' — 0 (44) 

($33/?al) / — ($44 + $45)0(o — ($45 +$55)060 “ 

(($13 + $36)060 ) ' ~ (($23 + $36)0>o ) ' = 0- (45) 

Upon integration we have: 

0oify) = F 0 (z) + z 
0oify) = Go(z)+z 
Pal(z) = 1+H 0 (z) 
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(46) 

(47) 

(48) 



and solving the imaginary part of the equilibrium equations: 


where 


H 0 (z) = 


z 

/ 


4>b\{z) = F 0 (z)-z 
ipbi{z) = Go{z)-z 
pb\{z) - l-Ho{z) 


2z + (Q13 + Q36) Fq + {Q23 + Q36) Go 
Q33 


dz. 


(49) 

(50) 

(51) 

(52) 


The displacement field in this case is: 


u x (x,y,z) = 

Ui ( x , y) Fo(z) + u 4 (x, y) z 

(53) 

Uy(x,y,z) = 

u 2 (x, y)G 0 (z)+u 5 (x,y)z 

(54) 

u z (x,y,z ) = 

m(x, y) + u G (x, y) Hq(z). 

(55) 

1.3 The model characterized by /3 2 


To find the mode of deformation for the model which satisfies the equilibrium 
equations up to the second power of we differentiate (23) to (25) twice with 

respect to (3 and let = 0. Upon integration, the following results 

are obtained: 

Mz) ■ 

— Fq(z) + z + Fi(z) 

(56) 

(z) ■ 

— Go(z) -h z + G\(z) 

(57) 

Pa 2 (z) 

= l + H 0 (z)-H 1 (z) 

(58) 

where 



aw = h 

( Mo Q 44 — No Q 45 u\ 

[ Q"Qx-c& Ho r 

(59) 
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R x dz 

\ Q44V55 — V45 ) 

(60) 
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(Q 13 + 2Qm + Q 23 )] _ 

(61) 

Q 33 J 
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and 


Mo{z ) = j + 2 Qi 6 + $66 — (Ql 3 + Q26) ** ^ Fo 
+ ^Qn + Ql 6 + Q26 + Qw - (Ql 3 + Q26 


-2z 


Qu + Q26I 


Q: 


dz. 


33 


( 62 ) 


N 0 (z) — J ^Ql2 + Ql6 + Q26 + Q66 — (Q23 + Qm) Fo 

+ (Q 22 + 2 Q 26 + Qm ~ (Q23 + Qm) ^ Gq 


-2z 


Q23 + Qoo 

Q33 


dz. 


Therefore the displacement field can be written in the form; 

u x (x,y,z ) = ui(x,y)F 0 (z)+u 4 (x,y)z + u 7 (x,y)Fi(z) 
u y (x,y,z ) = u 2 (x,y)G 0 (z)+u 5 (x,y)z + u 8 (x,y)Gi(z ) 
u z (x,y,z ) = u 3 (x,y) +u 6 (x,y) H 0 (z) + u 9 (x,y) Hi(z). 


(63) 


(64) 

(65) 

(66) 


This mode of deformation satisfies both the real and imaginary parts of the equi- 
librium equations up to the second power of (3. By continuing this process, the 
equilibrium equations can be satisfied to an arbitrary power of (3. 


2 The limiting case when — > 0 

One of the requirements of the hierarchic sequence of models is that each member 
converges to the same limit as the exact solution of the problem of elasticity as 
h — ► 0. We know from the evaluation of the laminated strip that the first two 
members of the hierarchy do not meet this requirement if some materials properties 
are not adjusted. We also learned that the adjustments are different for the first and 
for the second member of the hierarchy. The model which satisfies the equilibrium 
equations up to the second power of j3 is the first member of the hierarchic sequence 
of strip models which converges to the right limit as h — ► 0. 
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The exact solution minimizes the potential energy with respect to all functions 
Ui(x,y), t = l, 2, . . . for which the strain energy is finite. The limit for each model 
is obtained as follows: 

1 . Start with the expression of the potential energy: 



H-oo H-oo r+h/2 / 

J J J hj2 ^ x Gy ^~ x y ^ ^ xz ^ xz 

+Tyz~hz)dxdydz = j*™ J q(x,y)u x (x,y,h/2)dxdy 


(67) 


For a given plate model, compute the strain components and perform the integra- 
tion through the thickness [z direction) to obtain the material coefficients. These 
coefficients form the laminate material stiflhess matrix Eij. Rewrite the potential 
energy in terms of the Eij. 


U = ^ JJ^{u} T [E\ {u}dxdy- JJ^qu z (x,y, h/2)dxdy 

In the case of the /5° model, the potential energy expression is: 

dui dwi 


( 68 ) 




dy dy 




dx dx 
du 3 \ 2 
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dx J 
du 3 du 3 


dx J ^ dx dy 


+ Eg 


fdui\ 

U y) 

du 3 \ 
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fOU 3 \ 

U y) 


+ 


E10 ul+2EuUi ui + Euul +2 ^ 4 "^"^ + 

. „ du 2 dui , . „ du 2 du 2 , n du 3 du 3 . 0 r du 3 

2 ^^^ +2 ^ 6 ^^ +2E ' ,Ui -fc +2Eis ' u *~te +2Ei,u '~fy + 


2 E 10 u 2 dxdy — q(x,y)u,(x,y,h/2)dxdy 


(69) 


2 . Obtain the Euler equations by taking the variation of the potential energy with 
respect to each one of the field variables tt,(x,y), i — 1, 2 , — Apply Fourier 
transform to the Euler equations and construct the system of linear equations in 
the transformed field variables Ui(£,r)): 

[A] {C} = {R} (70) 

The matrix [A] depends on the material stiffness matrix Eij, and on the Fourier 
variables £, and 77 , and {i?} is the load vector obtained from the transformation of 
the potential of the external forces. 
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3. Compute values for Eij for different stacking sequences and solve the system of 
equations for U 3 for each one of them: 

D«, V ) U , = B(t, v ) Q (71) 

where 

D{(,V) = D l e + D,(r l + D,? + D t ? + D ! ,en + 

D s en‘ + Ih^ + D,r l t +-- (72) 


is the determinant of [A], and 

-£?(£> l) = -So + Bi + B 2 £r] + B 3 rj 2 H (73) 

is the determinant of [A] when the third row was replaced by the load vector {it!}. 
4. Perform the limit analysis as h — ► 0. First divide D{ by B 0 : 


A 
= A 

for i = 1 — ► 3, j — 1 — ► 3 

(74) 

Bq 

for i — 4 — ► 8, j — 1 — ► 5 

(75) 


Note that B 0 is the first non-zero term of f?(£, rf). Neglecting derivatives higher 
than fourth order we obtain the general equation : 

(<*1 f 2 + <*2 f V + 0:3 V 2 ) U 3 + 

(A £ 4 +Af 3 *?+A£V + A^ 3 +&V*) U 3 = Q. (76) 

This last equation has the following form after inverse Fourier transform: 


f &u 3 , cPu 3 
V* 1 dx 2 +a2 dx&y 


+ 0:3 


cPu 3 \ 

WJ 


+ 



+A 


(Pu 3 
dx 2 dy 2 


+ A 


&U3 

dx d ip 



= ?• 


(77) 


If the hierarchic plate model being evaluated converges to the proper limit, the 
coefficients a, should be zero. This is because all models should converge to the 
Kirchhoff model as h — ► 0 [3], which requires that a, = 0. If the a, are not zero, the 
model needs adjustment of the material properties. The /?, may require adjustment 
also, so that they have the same values as those of the models that converge to the 
same limit as the theory of elasticity as h — > 0. 

Following t his procedure, it has been found using Mathematica 1 that the model 
characterized by /3° is not a member of the hierarchy. However, making the trans- 
verse shear moduli constant through the thickness, the a, become zero. Further 
discussion of this model is deferred to next section. 


1 Mathematica: A system for doing mathematics by computer. Wolfrman Research Inc. 
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Figure 1: Model problems: Notation. 


3 Examples 

Consider a square plate of uniform thickness h and planar dimensions a and b = a, 
composed of perfectly bonded orthotropic layers, symmetrically distributed with 
respect to the middle plane, (Fig. 1). The number of layers, the stacking sequence 
(orientation of the material axes with respect to the global axes), and the boundary 
conditions are the variables used in several representative model problems. 

Uniform load is applied as a normal traction to the top and bottom surfaces of 
the plate. All layers in the laminate are of equal thickness, and are of a square sym- 
metric unidirectional fibrous composite material possessing the following stiffness 
properties, which simulate a high-modulus graphite/epoxy composite: 

E l = 25.0 x 10 6 psi Er = 1.0 x 10 6 psi 

G lt = 0.5 x 10 6 psi Grr = 0.2 x 10 6 psi 
vlt — vtt = 0.25 

where L indicates the direction parallel to the fibers, T is the transverse direction, 
and vlt is the Poisson ratio (i.e., vlt — —ctt/cll, where ejr, £ll are, respectively, 
the normal strains in the directions T and L). These material properties were 
selected from reference [1]. 
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When the L direction coincides with the x direction, we refer to it as the 0 = 0 ° 
orientation. For a three-plies laminate a designation 90 / 0/90 means that the central 
lamina is oriented with the L direction parallel to the global x-axis, and in the two 
outer layers L is at 90 ° with the global x-axis. 

In order to establish a reference solution which can be regarded as being suf- 
ficiently close to the exact solution of the problem of elasticity, the problem was 
solved using the finite element program MSC/PROBE 2 and an experimental pro- 
gram in which the algorithm described in this report is implemented. In the ref- 
erence solution obtained with MSC/PROBE each layer was discretized as a three- 
dimensional element with orthotropic material properties. The solution was ob- 
tained for p ranging from 1 through 8. in energy norm was below 1% at p — 8. The 
solution corresponding to p = 8 will be used as the basis for comparison. 

The solutions corresponding to the proposed hierarchic models were obtained 
using only one laminated plate element. The polynomial degree was varied from 1 
through 8 and the equilibrium equations were satisfied up to (3 1 . 

The model that satisfies the equilibrium equations up to the zeroeth power of 
/? was modified as indicated in the previous Section. The transverse shear moduli 

of each layer (Q 44 and Q 55 ) were made equal to the harmonic averages Qu and 

Q55, while Q45 was made equal to the average Q45. In the case of three layers for 
instance, the harmonic average of Q 44 is: 



The following changes were introduced for each layer: 

Qa = Qij - , i, j = 1 , 2, 6. ( 79 ) 

This modification is necessary to account for the plane stress constitutive equations 
for each layer used in the Kirchhoff and Reissner-Mindling type plate models [ 2 ], 
[ 6 ]. 

We will denote the modified model characterized by ( 3 ° with ( 3 ^. The plate 
deflection at a given location (x n , y„, z n ) was selected to evaluate the models. The 
normalized plate deflection is defined as: 

jj def 100 Et h? u z (x n , y n , 0) ( 80 ) 

qa 3 

where q is the applied traction, h is the thickness of the plate and u z (x n , y„, 0) is 
the vertical displacement of the middle plane of the plate at x = x n , y = y n . 

2 MSC/PROBE: User’s Manual, The MacNeal-Schwendler Corporation, 1600 S. Brentwood 

Blvd., Suite 840, St. Louis, Missouri 63144. 
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Normalized Central Displacement 
Three-plies Laminate (90/0/90) 


Uz 



a/h Ratio 


Figure 2: Orthotropic square plate. Four sides simply supported. 

3.1 Cross-ply laminate 

The results for a three-plies orthotropic (or cross-ply) simply supported square 
plate are shown in Fig. 2. For large a/h ratios both models yield similar results. 
A s a/h decreases, the /3° model underestimates the deflection while the /3 1 model 
is very close to the MSC/PROBE solution. 

Fig. 3 show the results for a three-plies (90/0/90) square plate with two opposite 
sides simply supported and the other two free. The central deflection of the plate is 
compared with the results of the laminated strip model described in reference [3] . 
Fig. 4 correspond to the same problem but with five layers (90/0/90/0/90). The 
agreement between the beam and plate models is excellent. 

The influence of the number of layers in the end deflection of a square plate with 
one end clamped and the other three free is shown in Fig. 5 for three different a/h 
ratios. In all cases the fibers in the outer layers were normal to the clamped edge 
of the plate. Also included in the figure are the results of the deflection computed 
using a simplified beam formula which is valid for the case a/h — *• oo. According 
to reference [4] the end deflection of a cantilever beam of length a and thickness h 
with uniform load q is: 
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Figure 3: Orthotropic square plate. Two sides simply supported. 


Normalized Central Displacement 
Five-plies Laminate (90/0/90/0/90) 
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Normalized End Displacement 
Influence of Number of Layers 


Uz 


3 30 300 
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Outer fibers normal to damped edge 

Figure 5: Orthotropic square plate. One side damped. 
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where 


+h/2 

D n - f Qn z 2 dz 
-h/ 2 


(82) 


The deflection computed by the use of (81) is identified as ‘Beam (a/h — ► oo)\ The 
results indicate that when the number of layers increases the bending stiffness of 
the plate decreases to an asymptotic value. In the limit when we have an infinite 
number of alternating plies, the laminate will become quasi-homogeneous. The 
property of the laminate will be square symmetric but not homogeneous [5] . 


3.2 Angle-ply laminate 

The results for a three- and four-plies simply supported square angle-ply laminated 
plate are shown in Fig. 6 and in Fig. 7 respectively. In this case the stacking 
sequence is such that all layers are oriented at either 45° or —45°. For the three- 
plies plate, the agreement between the MSC/PROBE solution and the hierarchic 
models is similar to the case of cross-ply laminates. For the four-plies plate no 
MSC/PROBE results are still available. 
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Normalized Central Displacement 
Three-plies Laminate (-45/45/-45) 


Uz 



Figure 6: Angle-ply square plate. Four sides simply supported. 


Normalized Central Displacement 
Four-plies Laminate (-45/45M5/-45) 


Uz 



Figure 7: Angle-ply square plate. Four sides simply supported. 
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Normalized Central Displacement 


Uz-2 Uz-4 



Central Layer Angle (deg) 


Results for a/b»10 
7 *fr 2-sides SS (Beta Om) 
g-sides $$ (Beta 1)' ; 
4-sides SS (Beta Om) 
•B" . 4-sides $S (Beta 1) 


Figure 8: Angle-ply square plate. Influence of fiber orientation. 

The influence of fiber orientation in the central deflection of a three-plies square 
plate with two opposite sides simply supported (2-sides SS) and Four sides simply 
supported (4-sides SS) is shown in Fig. 8. In this case the a/h ratio was kept 
constant at a/h = 10 an the orientation of the fibers in the central layer was varied 
between 0 and 90°. The fibers in the outer layers were always at 90° with the ones 
of the central layer. The results for the /3° and j3 l models are included for each 
boundary" condition. 


4 Summary and conclusions 

1. Hierarchic models for mid-plane symmetric laminated plates have been de- 
veloped based on a single parameter /3. The powers of the parameter (3 rep- 
resenting the degree to which the equilibrium equations of three-dimensional 
elasticity are satisfied, have been used to identify members of the hierarchic 
sequence. 

2. The model characterized by /3° is the Reissner-Mindlin model, generalized for 
laminated composites, when the modified material properties are used (also 
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known as first order shear deformation model). In the special case, when 
the shear modulus is independent of z, the hierarchic model is the Reissner- 
Mindlin model. The shear correction factor can be assigned arbitrarily since 
the requirements set for hierarchic models are satisfied independently of the 
shear correction factor. 

3. The advantage of a single parameter sequence of models is that the number of 
fields added per level is always the same. Three fields are added per increment 
of power of /?, regardless of the number of layers. 

4. Good correlation between the proposed hierarchic sequence and a three- 
dimensional solution has been found for the problems investigated. More 
work is under way to compare not only displacements but stress distributions 
as well. 
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